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Abstract 

Background: The protein p53 has a well established role in protecting genomic integrity in human cells. When 
DNA is damaged p53 induces the cell cycle arrest to prevent the transmission of the damage to cell progeny, triggers 
the production of proteins for DNA repair and ultimately calls for apoptosis. In particular, the p53-Mdm2 feedback 
loop seems to be the key circuit in this response of cells to damage. For many years, based on measurements 
over populations of cells it was believed that the p53-Mdm2 feedback loop was the responsible for the existence 
of damped oscillations in the levels of p53 and Mdm2 after DNA damage. However, recent measurements in 
individual human cells have shown that p53 and its regulator Mdm2 develop sustained oscillations over long 
periods of time even in the absence of stress. These results have attracted a lot of interest, first because they 
open a new experimental framework to study the p53 and its interactions and second because they challenge 
years of mathematical models with new and accurate data on single cells. Inspired by these experiments standard 
models of the p53-Mdm2 circuit were modified introducing ad-hoc some biologically motivated noise that becomes 
responsible for the stability of the oscillations. Here, we follow an alternative approach proposing that the noise 
that stabilizes the fluctuations is the intrinsic noise due to the finite nature of the populations of p53 and Mdm2 
in a single cell. 

Results: We study three stochastic models of the p53-Mdm2 circuit. The models capture the response of the p53- 
Mdm2 circuit in its basal state, in the presence of DNA damage, and under oncogenic signals. They are studied 
using Gillespie's simulations, mean field methods and analytical predictions within the context of the Linear Noise 
Approximation. For the first two models our results compare quantitatively well with existing experimental data 
in single cells. The study of the response of the p53-Mdm2 circuit under oncogenic signals is process for which 
we do not have single cell measurements, but can be modeled waiting for experimental verifiable results. While 
we can not discard that other sources of noise in the cells may also be important, our results strongly support the 
relevance of the intrinsic noise in these systems. 

Conclusions: We suggest that the noise induced by the finite size of the populations is responsible for the existence 
of sustained oscillations in the response of the p53-Mdm2 circuit. This noise alone can explain most of the 
experimental results obtained studying the dynamics of the p53-Mdm2 circuit in individual cells. 
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Background 

The p53 protein has attracted special attention dur- 
ing the last twenty years because there is strong ev- 
idence of its malfunction in most human cancers [T] . 
Like most proteins, the p53 may be inactivated di- 
rectly by mutations in its own gene or as a result of 
the interactions with other genes that transmit infor- 
mation to or from the p53. Unfortunately, the many 
pathways that control its activation and the com- 
parable large number of functions, some apparently 
contradictory, still makes the complete understand- 
ing of the role played by this protein an unsolved 
problem [TJ[2]. 

In general, p53 works as a transcription factor 
that positively or negatively regulates the expres- 
sion of several, and very different genes. The p53 
gene was first identified as a tumor suppressor gene 
already many years ago, but we have now evidence 
that it has a role in the regulation of glycolysis [3] , 
the repair of oncogenic response pQ, the regulation 
of metabolism [5] and many others [5] . 

It is now well accepted that one of the main 
mechanisms regulating the expression of the p53 is 
its interaction with the protein Mdm2. The p53 and 
the Mdm2 form a feedback loop very common in 
many biological systems. In normal conditions the 
activation of the p53 activates the Mdm2 that in 
turn suppresses the p53 |6]. Then, if the cell is in 
the presence of stress, like DNA damage, p53 and/or 
Mdm2 are phosphorilated and its interaction is re- 
duced increasing the level of p53 in the cell [6]. 

An important breakthrough in the study of this 
system was achieved a few years ago when re- 
searchers tracked, cell by cell the concentration of 
p53 and Mdm2 after irradiation with gamma rays 
[7]. More recent experiments of the same groups 
were fundamental to clarify at least the following is- 
sues [7H9]: 

• The dynamics of the p53 is dephased from the 
dynamics of the Mdm2 and both are charac- 
terized by a fixed frequency and variable am- 
plitudes of oscillations within a single cell. 

• The dynamics and the quantity of the protein 
p53 is similar before and after irradiation. This 
open the question of how the p53 leads to cycle 
arrest or apoptosis when the cell is damaged. 

• The levels of p53 are uncoupled from its tran- 
scriptional activity. 



From the mathematical point of view these re- 
sults inspired the development of new models. Most 
of them were built on previous ideas about the in- 
teractions within the p53-Mdm2 circuit but adding, 
ad-hoc, some biologically motivated noise that will 
be responsible for the stability of the oscillations 
[ZHiniHI]. Others, introduced a delay mechanism 
between the p53 activation and the appearance of 
the Mdm2 [HE3]. 

On the other hand, experiments suggest that the 
size of the population of p53 in the basal and ac- 
tivated states are rather small, 10 3 — 10 4 elements 
(see for example [13]) suggesting than mean field ap- 
proaches may not be appropriate. Here, we follow an 
alternative approach and propose that the finite size 
of the populations involved is the source of an intrin- 
sic noise that, in a single cell, stabilizes the oscilla- 
tions. This kind of approach was recently suggested 
in ref. p3] but here we push it further. First studying 
not only the isolated p53-Mdm2 module but also the 
module in the presence of different external stresses: 
DNA damage and oncogenic response. Second, us- 
ing analytical arguments we are able to clarify the 
role of the different parameters of the models in the 
stability of the fluctuations. Moreover, we show that 
this mechanism consistently reproduces some of the 
experimentally available data. 

The p53-Mdm2 model and Mathematical Back- 
ground 

Probably the simplest version, and the starting point 
of any development of the p53-Mdm2 module is the 
loop presented in Figure 1. This is a feedback loop 
in which the expression of p53 activates the pres- 
ence of the Mdm2 that in turn regulates the activ- 
ity of p53. This is a well studied model |7|, whose 
mean field solution is characterized by damped oscil- 
lations. Sometimes, it is convenient to introduce an 
intermediary for the interaction between the p53 and 
the Mdm2, usually a messenger mRNA that acts as 
the transcription factor of the Mdm2. This interme- 
diary is the responsible for a dephasing between the 
oscillations of the p53 and Mdm2 proteins [15] but 
for simplicity we will assume in the rest of the work 
that it is in equilibrium, keeping, as our basic block 
for the p53-Mdm2 module, the scheme of Figure 1. 

We formalize the model starting from a stochas- 
tic approach considering the number of components 
of each species, p53 (71^53) and Mdm2 (jiMdrra) as 
the relevant variables of the problem. 
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From this point of view one faces the problem of 
solving the following Master equation: 

(i) 

where P{n, t) is the probability that the system is 
in state ft = {n P 53,nMdm2) at time t and W n ^~, are 
transition rates from one state to the other. The so- 
lution of eq. (TTJ) is not an easy task, but, for very 
simple transition rates and one most usually use nu- 
merical methods, or resort to approximations. In 
this work we try both approaches and compare the 
results. 

The numerical solution is obtained using the 
Gillespie's algorithm [16]. As other Monte Carlo 
techniques this algorithm guarantees an exact solu- 
tion, but it is time consuming and is not a good tool 
to explore the parameter space. To understand the 
role of the parameters of the model in the sustain- 
ability of the fluctuations, it is important to have at 
least approximate solutions to the problem. Here, 
we use the Linear Noise Approximation (LNA) [T7] , 
that was first proposed in the context of chemical 
kinetics and have gain considerable attention in the 
last few years for modeling intracellular processes, 
[HUH], but also in more general contexts [201121] . 
The mathematical details of this approximation can 
be found in [17] . and we summarize them in section 
Methods. However, to simplify the reading, and to 
keep the consistency of the manuscript we present 
here its main results. 

One of the difficulties to solve equation (fTJ is the 
discrete character of rii. To deal with this, van Kam- 
pen |17j proposed the following approximation: 

rii = Qxi + VTicti (2) 

i.e., one approximates the number of molecules of 
specie i by some concentration, plus some fluctua- 
tions proportional to the square root of the system 
size f2. Substituting @ into JT]) and developing the 
equation in order of ^ one obtains, at first order a 
set of differential equations for xi , usually called the 
mean field or deterministic approximation: 

3 

where Sij is the stoichiometric matrix of the interac- 
tions (see the section Methods for a detail deriva- 



tion). At second order one can prove that the fluctu- 
ations cti are Gaussian variables. However, to com- 
pare them with the experiments it is more useful to 
describe the evolution of these variables in terms of 
Langevin equations. 

a N 

— = 2^A ik a k + rn{T) (4) 
fc=i 

where A ik = J2jLi S ik^f- and the ^ has zero 
variance and correlation function (i]i( T ) r lj( T ')) = 
DikS(r - t) with D ik = Y,jLi SijSkjWj(x) (see 
Methods). 

Equation Q is a linear Langevin equation that 
can be easily solved in the Fourier representation 

Q.W =^Kk lf lk(w) (5) 

k 

where = iwSik — A ik . From a^w) one can get 
the power spectrum of the fluctuations as: 

P t {w) = (I5.HI 2 ) =EE^ 1 *^ 1 ^ ( fi ) 

3 k 

which, when studying fluctuations in single cells, is 
the quantity usually reported experimentally 0[8]. 

Results and Discussion 

In this section we present the main results of our 
work. It is divided in three parts, in each one of 
them the p53-Mdm2 module is studied either as an 
isolated structure or considering its interactions with 
different external elements. The organization of each 
subsection is the same, we first motivate the model 
from the biological point of view and then present its 
stochastic formulation. We will show the results of 
simulations using Gillespie's algorithm and compare 
them with the prediction of the corresponding mean 
field solution ([3]). Then, we study the fluctuation 
spectra obtained from the simulations and compare 
them with the predictions of the Linear Noise Ap- 
proximation (LNA). Finally when it is possible we 
discuss the connection with the experimental results. 

Basal response 

In the absence of any external signals the p53, is 
activated only by random DNA damage. This is a 
common event during the cell cycle life. Spontaneous 
hydrolysis, collapsing replication forks and oxidative 
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stress [2"2l2"3"] , but also metabolic stress may produce 
p53 activation [5] . Within this context we study the 
basal response using the model in Figure 1. 

In stochastic terms the model is well described 
by the following set of expressions: 



p53 



Wi 



p53 

Wi 



p53 

2>53 
Mdm2 



Mdm2 



2p53 



II'h 



4> 



p53 + Mdm2 —4 Mdm2 



(7) 
(8) 
(9) 
(10) 

(11) 



where equations ([7]) and ([8]) reflect the activation of 
Mdm2 and the self-activation of p53. Equations (J9j) 
and ([TU)) the decay of both species and equation (fTT|) 
reflects the regulation of p53 by Mdm2. The specie 
$ is included to allow the free variation of the quan- 
tities of p53 and Mdm2 during the simulation. The 
different Wi are transition probabilities proportional 
to the size of the system [17]. To fix concepts we 
fixed the simplest possible transition probabilities, 
in this case: 



Wi = kxflXp^ 


(12) 


W 2 = k 2 flx p53 


(13) 


W 3 = k 3 Qx p53 


(14) 


Wi = k i riXMdm2 


(15) 


= k 5 £lx p 53XMdm2 


(16) 



One should notice that the transitions (Wi,2), i-e, 
activation of Mdm2, and self- activation of p53 are 
proportional to the quantity of p53 in the system. 
The degradation of Mdm2, (W4) is proportional to 
the amount of Mdm2. The p53 is degraded following 
two mechanisms, a normal degradation (W 3 ) pro- 
portional to the quantity of p53 in the cell, and a 
degradation activated by Mdm2 (W5), that is pro- 
portional to the quantities of p53 and of Mdm2. 

Following ([3]) the mean-field equations of this 
model read: 



dx 



p53 



df 



= k 2 Xp 53 - k 3 X p53 - k 5 X p53 XMdm2 (17) 
dXMdm2 



dt 



= kiX p53 - kiX M dm2 (18) 



The fixed points of these equations are: (0, 0) and 
(flit' If) wnere to simplify notation we considered 



k' 2 = k,2 — k 3 and relabeled it as k 2 . Then, the non 
trivial solution has biological sense only if the rate of 
activation of p53 is larger than its degradation rate 
(fc 2 >0). 

The stability of these fixed points is defined by 
the following Lyapunov exponents: 



(0,0) A 



k 2 

—k 4 



(19) 



and: 



p53' X Mdm2 



)->A 



i(-fe 4 + y /k'j - Ak 2 k A ) 
i(-fe 4 - \/k\ - Ak 2 k 4 ) 



(20) 

These results guarantee that if the non-trivial so- 
lution exists ( &4 > and k 2 > ), then it is stable. 
Moreover, if k 4 < Ak 2 , i.e. if the activation of p53 is 
large enough, it displays damped oscillations. This is 
the region of parameters relevant from the biological 
point of view. 

In Figure 2 we show the comparison between the 
numerical solution of these equations with the results 
of the Gillespie's algorithm for the same values of the 
parameters. Note that, while the deterministic solu- 
tion is over-damped, the stochastic solution displays 
persistent oscillations around the fixed point. 

To characterize these sustained oscillations we 
computed the power spectrum of the results of 
the Gillespie's simulations over 200 realizations and 
made their geometrical average. The resulting power 
spectrum of the oscillations of the p53 is presented 
with points in Figure 3 where we also show the an- 
alytical results predicted by the LNA (see equation 
(J6j) and section Methods) . This plot compares very 
well with Figure 6 in [8] where the authors studied 
the basal response of the system and presented ex- 
perimental results for the Fourier transform of the 
fluctuations. 

The existence of a maximum in the power spec- 
trum of the fluctuation in Figure 3 is associated to 
the presence of a characteristic frequency for these 
fluctuations. In an infinite system, these oscillations 
die out, but the stochasticity produced by the finite 
size of the system acts as a white noise, see equation 
(j4{, that interacts with the proper frequency of the 
system producing a resonance-like effect. This char- 
acteristic frequency is connected with the presence 
of an imaginary part in the eigenvalues of the system 
of differential equations. 

This is made evident by a simple inspection of 
the LNA expression for the power spectrum of the 
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model (see Methods). There, one can easily check 
that the denominator has a minimum, very close to 
the expected frequency of the damped oscillations 
estimated in the mean field model. Below, we show 
that this phenomenon repeats in the other models, 
although in a more sophisticated form. 



Response in the presence of stress 

The ATM is a kinase protein that is currently in- 
volved in the detection of the damage in the DNA. 
When the ATM detects the damage it induces the 
phosphorilation (activation) of p53. We model the 
presence of this kind of external stress incorporat- 
ing to the module of the basal state the new protein 
(ATM) and its interaction with the p53 [21E]. 

In Figure 4 we present a sketch of this model. 
From the technical point of view, we have now a 
problem where to the dynamics governing the basal 
state we must add the following new transitions: 



Hi 



ATM ATM + p53 



ATM + p53 p53 
ATM ^ $ 



ATM 



2ATM& 



(21) 
(22) 
(23) 
(24) 



These transitions reflect the activation of p53 by 
the ATM, (We), the regulation of the ATM by the 
p53, (W7), the natural degradation of the ATM, 
(Ws), and a self-reinforcement field for the ATM, 
(Wg) that mimics the persistence of the damage. As 
in the previous section we assumed that they take 
the following simple form: 



W 6 = k 6 £lxATM 
W 7 = k 7 ttx p53 x A TM 
W s = k 8 ttx ATM 

Wg = kgttx AT M 



(25) 
(26) 
(27) 
(28) 



With these definitions of the transition rules, we 
may write the mean field equations of the model: 



dt 



dx 1 



dt 



k2X p53 - k 3 x p53 - 
ksXpZzXMdml + kgXATM 
dXM J™ 2 = k lXp 53 ~ UxMdral 
= kgXATM - k S X A TM ~ k 7 X A TMX p $ 3 



(29) 



These equations are similar to the ones in (TI"5|) . 
but now the role of ATM on the activation of p53 
results in the addition of a new term to the first 
equation. Moreover, ATM is also a new dynamical 
variable whose evolution must be taken into account. 

This system of equations has three different 
fixed points. They are, (0,0,0), (j^, fj, 0) and 

(&> & " **)) where for simplicity we 

re-labelled k' 2 — k 2 — k 3 , and k' 9 = kg — ks, as ki and 
kg respectively. 

The Jacobian matrix necessary to study their 
stability takes the general form: 



k 2 k 5 x* Mdm2 
ki 

—k 7 x* ATM 



that evaluated 

( kc) kyki ko l kcjk\k 

V k-Y ' ^7^4 ' kjks V 

lowing secular equation: 



in the fixed 



^5^p53 


k 7 


— &4 











point 


of 



interest 



£2)) transforms into the fol- 



A 3 + A 2 (c-fc 4 ) + A( 



kgk x k 5 



k- 



where c = k 2 - 



k^c — kgc) — kgk^c = 
(30) 



k-jki 



This secular equation must be studied numeri- 
cally and its solution defines the stability of the fixed 
points. We have found that depending on the param- 
eters, three different regions are well defined. In Re- 
gion I the solution {x* 53 , x* Mdm2 , x* ATM } may exist, 
is stable and the system reaches the fixed point with 
damped oscillations. In this region also the solution 
{x* 53 ,x* Mdm2 ,0} is stable if k 2 > 0. Depending on 
the initial condition, the system is attracted to one 
or another fixed point. 

In Region II the solution {x* p53 ,x* Mdm2 ,x* ATM } 
has biological sense and is stable, but the fixed 
point is reached without damped oscillations. In 
this region the solution {x* 53 , x* Mdm2 , 0} does not 
exist. Finally, in Region III only the solution 

{ x p53' x *Mdm2> 0} is stable. 

On the basis of the results obtained study- 
ing our previous model we concentrate our atten- 
tion in Region I, where the mean field fixed point 
{Xp 53 ,x* Mdm2l x* ATM } is reached through damped 
oscillations. In Figure 5 we show the numer- 
ical solution of equation (|30[) and the results of 
the Gillespie's algorithm using the parameters k — 
(0.99, 1,0.44, 0.69, 0.85, 0.5, 0.5, 0.1, 0.^/i" 1 . As in 
the previous section, we see that the stochastic sim- 
ulation shows sustained oscillations in time. 
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To characterize these oscillations we use again 
the power spectrum of the Gillespie's results. In Fig- 
ure 6 we show with points this power spectrum and 
with a continuous line the predictions from the LNA. 
Note again the good coincidence, but one can also 
compare them with Figure 3 in [7]. Our model re- 
produces qualitatively and quantitatively the experi- 
mental observations: the existence of a characteristic 
frequency and a divergence close to zero frequency. 

The peak in the power spectrum still reflects the 
resonance-like mechanism between the p53-Mdm2 
feed-back loop and the intrinsic noise. On the other 
hand, the divergency at low frequencies results from 
the presence of the new interaction and can be in- 
terpret looking to the analytical form of the power 
spectrum derived in the LNA at low frequencies. It 
has the form 

P p53 (^^0)- ,,, 2fc5fc I (31) 

KlK 6 Ki 3 - k 2 k 5 k 9 

and guarantees large values of the power spectrum 
if the denominator is small enough. 



and defining the transition rates as: 

W 6 = k 6 Qx p i4,XMdm2 (36) 

W 7 = kjQxp^Xpn (37) 

W 8 = k & rtx pU (38) 

W 9 = k 9 Qx pl4 (39) 

we obtain the following system of equations for the 
dynamic of the system: 

-^p = k 2 X p53 - k 3 X p53 - k 5 X p53 X M drn2 (40) 
+ k 6 X A TM 

dXA ^ m2 = kiX p53 - kiXMdm.2 - k 6 X p i4X M dm2 

= k 9 x p u - k 8 x p u - k 7 Xpux p53 

which again depending on the parameters may 
show three fixed points: (0,0,0), (f^, ia.,0) and 

(fe>!f< khrS^tr ~ W) where to simplify the 
notation we used k 2 = k 2 — fc 3 and k 9 = k 9 — k-j. 

The Jacobian matrix near these fixed points 
takes the form: 



Oncogenic response 

Another important pathway for the activation of 
p53, comes from the activation of the plA ARF pro- 
tein. In the presence of oncogenes this protein is 
activated and joins the Mdm2 accelerating its degra- 
dation [24 . This process leads to the increase of the 
p53 level in the cell. 

As far as we know there is no experimental evi- 
dence for the single cell dynamics of p53 and Mdm2 
in the presence of oncogenes, therefore the parame- 
ters used here have no experimental bias. Neverthe- 
less we will show that independently of the values 
used the power spectrum of the fluctuations is qual- 
itatively different from the one obtained when the 
p53 interacts directly with ATM. 

The model of oncogenic response is sketched in 
Figure 7. To the standard interactions present in 
the basal model we add the regulation of pl4 ARF to 
Mdm2 and a regulation of pl4 ARF by p53. In terms 
of stochastic transitions we have now: 



Mdm2+pU ARF ^ P U ARF 
P U ARF +p53^\p53 

pU ARF $ 
pU ARF 2pU ARF^ 



(32) 

(33) 
(34) 
(35) 



J = 



■ k§x* Mdrn2 

ki 
-k 2 x* 14 



-k^x 



p53 
A,'4 





> Jj Mdm2 
kg 



(41) 

that evaluated in the fixed point {x* 53 ,x^ ldm2 ,x* 14 } 
leads to the following secular equation: 



A 3 + A 2 fc 4 + A( 



k 9 hk 5 kgk 5 ki 



k 2 k 9 k 4 = (42) 



which defines a phase diagram qualitatively simi- 
lar to the one discussed in the preceding section. 
Three regions, but only in one of them, the so- 
lution {a;* 53 , x* Mdm2 , x* 14 } displays damped oscilla- 
tions. Our previous results suggest that this is the 
relevant region for our purposes. 

Then , using as parameters k = 
(1.0,1.0,0.44,0.69,0.85,0.5,0.5,0.4) we obtain 
damped oscillations in the mean field solution of 
the problem and sustained oscillations through the 
Gillespie's algorithm. The results are shown in Fig- 
ure 8. 

In Figure 9 we show the power spectrum ob- 
tained using the LNA and the one obtained from the 
simulations. Also for this model, the LNA approx- 
imation describes correctly the results of the simu- 
lation. As before a clear peak representative of the 
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feed-back loop between the p53 and Mdm2 is ob- 
served, however, there is no evidence of a divergence 
at small frequencies. 

This absence of the divergence is independent of 
the parameters chosen. It can be understood look- 
ing again to the structure of the power spectrum 
predicted by the LNA at small frequencies. For this 
model: 



Pp53(w -> 0) 



2fc 5 (fci + fc 2 ) 
klk 6 



(43) 



that indicates that only when k\ or fc 6 are zero, i.e, 
only when the regulation of the Mdm2 by the p53 
or the interaction between the Mdm2 and plA A RF 
disappears we have a diverging power spectrum at 
small frequencies. But in this case our model of in- 
teractions breaks down. 



Conclusions 

We studied, using analytical tools (LNA) and nu- 
merical simulations (Gillespie's algorithm) the p53- 
Mdm2 module in its basal state, in the presence of 
DNA damage and oncogenic signals. The module 
is studied as a stochastic system where fluctuations 
are relevant. Our results compare quantitatively well 
with the experimental data available. This suggests 
that the intrinsic noise due to the finite size of the 
populations of p53 and Mdm2 may be responsible 
of the sustained oscillations in the response of the 
p53-Mdm2 circuit. 



Methods 

Linear Noise Approximation 

The quantities characterizing the number of parti- 
cles in our problem are discrete. The deviation of 
this number from the expected mean value produces 
fluctuations usually called intrinsic noise because it 
is not caused by any external source. When the sys- 
tem is small enough, this intrinsic noise is not neg- 
ligible, and a mean field approximation can not be 
taken. 

In these cases, the correct way to describe the 
system is to solve the Master Equation that governs 
it: 



dP(n,t) 
Ft 



^(P(n',t)W n ^ n P(ft,t)W n ^,). 

(44) 



This is an equation for the change in time of the 
probability P(n,t) of being in the state n at the 
time t, and it rests upon the Markovian hypothesis. 
W~ , -, are transition rates from states n to n'. 

An important difficulty in solving the Master 
equation arises from the discreteness of the variables 
involved (ft). Its solution is obtained by analytical 
methods in very lucky situations. An approximation 
scheme that allows to get some information from it, 
is the linear noise approximation (LNA). 

This approximation assumes the possibility of 
studying the number rij of integrants of the i th pop- 
ulation as a main part Qxi that represents the mean 
value of the population size, plus a small deviation 
VUai that represents the fluctuations around the 
mean value, and that is proportional to the square 
root of the system size il. In this way: 



(45) 



As every process that takes place in the system 
causes variations of one or more integrands in the 
populations of the species involved, it is useful to 
write the master equation in terms of these discrete 
fluctuations. For this purpose let us define the step 
operator as: 

£?/(..., n h ... ) = /(..., m + k,...). (46) 
In terms of this operator the master equation reads: 



dP(n, t) 
~dt 



M 



N 



Wj(n,tl)P(n,t), 
(47) 

where Wj is the microscopic reaction rate for the j 
process and SV,- is the element of the stoichiom- 
etry matrix S, the magnitude of the change in the 
i th population when the j th process takes place. The 
total number of different populations that form the 
system is N, and M is the total number of different 
kinds of processes that can take place in the system. 

In terms of the fluctuations, which are continu- 
ous variables, the step operator takes the following 
form, more suitable for an analytic treatment 



fc _LA fc2 92 



+ ... 



(48) 



in this way, we can work with a continuous opera- 
tor, and still not lose the stochastic features of these 
processes. With it 
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JV 
i=l 

JV JV 



JV 



H 5jj 5fc j - „ 4 



i=l fc=l 



9a; 9a/, 



(49) 



Now we are going to replace the probability 
P(n, t) by the probability II(a, t) of being ^/Tid away 
from the mean field prediction 



p(n, t) = p{Qx + voa) = n(a, *) (50) 

We know that: 



dU(d, t) dP(n, t) n A 9P(n, i) 9a;, 

'- S 2 >^ — — and 



9f 



dt 



dm dt 

dP(n,t) _ 1 dU(a,t) 
9n, t/Q, dcti 

with this results we can obtain the relation 



JV 



dU(a, t) i ^ dU(a, t) dxi _ 8P(n, t) 



dt 



9ai 9t dt 



(51) 



If now we expand the transition probabilities per 
unit time around x 



JV 



dxi 



aj + o(fi x ) 



(52) 



Substituting (|4^|) and ([5"2"j) in the master equation 
(|4"T)l and using (fST))) . after grouping the coefficients 
of the similar powers of f2 we find: 



If now we compare the coefficients of the similar 
powers of £1 in equations (|53|) and ([5T|). where the 



unit of time has been changed from t to r 
find for the leading order: 



77 



M 

E 



^I^-(x). 



wo 



(53) 



This corresponds to the deterministic rate equa- 
tions that are often used to describe the dynamics of 
such systems at a mean-field level. It has been dis- 
cussed already in the subsection referring to mean 
field equations and stability analysis. 

For the coefficients of f/! -1 the relation reads: 



gn(q,r) 

97 



d[IL(a,t)a k 

Aik 



i,fc=l 



1 ^ 



Da 



9a, 
d 2 U(a,t) 



(54) 



^ 1 doodah 

i,k=l 



where 



A I 



Aik — ^ ' S. 

3 = 1 



dw 3 (x) 



dx k 



(55) 



is the jacobian matrix of the system, and 

M 

Dik = E S ij S kjWj(x). 
3=1 



(56) 



Equation (|54l) is a Fokker-Plank equation for the 
probability for the system to have the deviation y/fid 
around the mean-field prediction. This equation 
can be written in a completely equivalent formula- 
tion more benevolent to investigation using Fourier 
transforms. The problem can then be formulated 
as the set of stochastic differential equations of the 
Langevin type: 



dt 



= n (E^)E^^^- 

9[±I(a,t)a fc ] dWj(x) 



3=1 i,k=l 
M JV 



-n- 1 ( E E * 

, 3=1 i,fe=l 



dxk 



M JV 



I EE^ H 7-' /; 

3=1 i=l 



^]+o(fl-») 



77 = E 

fc=i 



(57) 



where %(r) is a Gaussian noise with mean and 
correlation given by (^(r)^^')) = D lk 8(T - r'). 

This a linear Langevin equation that can be eas- 
ily solved in the Fourier representation 



(58) 



s 



where <$>ik = iwSik — A^. From ai(w) one can get 
the the power spectrum of the fluctuations as: 

p,(«0 = (la.Mi 2 ) = EE^ 1 **^ (59) 



Power spectrum in the basal state 

Now we are going to analyze the general expression 
for the power spectrum (|S"9")l for the basal state. The 
analytics for the other two models develops following 
the same steps, with a more cumbersome Algebra 

N N 

Pi( U ) cx EE^H^^hH (60) 

3=lfc=l 



the denominator becomes small, and the power spec- 
trum has a large peak centered on this frequency. 
The condition dPi(uj)/duj — gives: 



2(w 2 - det(A))2oj + 2uj{Tr{A)) 2 = 



(64) 

2 _ 



Then the condition (|64[) simply becomes: oj 
det(A) - ( Tr ^)) 2 ; it implies that: 2det(A) > 
(Tr(A)) 2 , this condition in terms of the stability ma- 
trix implies that the eigenvalues of A are complex. 
It means that if the power spectrum has an extreme 
the stability matrix has complex eigenvalues. 

Replacing in the expression for the power spec- 
trum the eigenvalues for the model of the system 

in basal conditions: Xi 



-fc 2 

2 



fc 2 fc 4 



k 2 

4 



and 



a 2 = =!j* - iy / fc 4 fc 2 - \ 



Since $ = icuSik—Aik and since A and D are inde- 
pendent of oj, Pi{uj) is a fraction, the numerator and 
the denominator have both a polynomial structure 
of order 2N . The explicit form of the denominator 
is: \det$(oj)\ 2 . 

Then, we expect that the enhancement of the 
fluctuations would be in the values of oj that mini- 
mize the denominator. Now we are going to analyze 
the general expression for the denominator of the 
power spectrum, and show explicit expressions for 
the basal response. 

The matrix A in the expression of $ is the ja- 
cobian matrix which can be written in terms of the 
eingenvalues, then we can write: 



det^ (oj) — det 



Pi(u) cc 



iui — Ai 



-lu 2 - iuj(\i + A 2 ) + AiA 2 



(61) 



Pi(u) oc 



(—oj 2 — iui(\i - 


-A 2 ) -+ 


■ AiA 2 ) 


1 






(—a; 2 + iui(Xi H 


-A 2 ) 4 


- Ai A 2 ) 


1 






( UJ 2 -\ 1 \ 2 ) 2 + 


oj 2 (\i 


+ A 2 ) 2 


1 







(62) 



(63) 



(oj 2 - det(A)) 2 + oj 2 (Tr(A)) 2 



This form for the power spectrum shows clearly 
the existence of a resonance: for a specific value of oj 2 



P 1 (oj) 



Pi M « 



1 



cj 4 + cj 2 (fc 2 -2fc 4 fc 2 ) + /c|fc 2 
4oj 3 + 2oj(k\ 



2k A k 



4«2j 



(oj i + oj 2 (k 2 -2k i k 2 ) + k 2 k 2 ) 2 



(65) 



= (66) 



Pi(uj) has an extreme in oj = <J k A k 2 — ^ 
As we can see, the imaginary part of the eigen- 
values (7m[Ai j2 ]) and the value of frequency where 
Pi(oj) has an extreme are similar. Now we are going 
to write it in a more clearly way: 



Im[\x t 2\ = \ k 2 k 



k 2 



4fc 2 fc 4 



8fc 2 fc4 



(67) 

The above Taylor expansion may be a good ap- 
proximation if we take into account that for the cor- 

responding values of the parameters 4fc2 " fc4 < 
0.5. 

In the same way: 



2&2&4 



< 



k/\k 2 — = \l kokA \ 1 1 — 

2 V V 2fc 2 fc 4 



ft/4 

4fc 2 fc 4 ' 



\A2Ml 



Summarizing, the power spectrum of the fluctu- 
ations has a maximum at a frequency close to the 
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imaginary part of the eigenvalue that characterizes 
the frequency of the damped oscillations in the de- 
terministic regime. 
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Figures 

Figure 1 - p53-Mdm2 basic module 




Sketch of the p53-Mdm2 feedback loop in its basal state. The protein p53 activates Mdm2 and Mdm2 
suppresses p53. 
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Figure 2 - Oscillations of the p53 in the basal state 
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time (hours) 

Typical run showing the oscillations of the p53 in the basal state. The bold curve indicates the mean field 
solution and the thin one are results from the Gillespie's simulation. Parameters: k = (0.99, 1, 0.44, 0.49, 1.05)/i 
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Figure 3 - Power spectrum of the fluctuations of the p53 in the basal state 
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frequency f (1/h) 

Average power spectrum of the fluctuations of p53 in its basal state. The points represent the power 
spectrum obtained after averaging over 1000 realizations of the Gillespie's algorithm. The continue line is 
the prediction from the LNA. Note the existence of a characteristic frequency of oscillations. Parameters: 
k = (0.99, 1, 0.44, 0.49, lM)^ 1 . 
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Figure 4 - p53-Mdm2 model in the presence of DNA damage 




Sketch of the p53-Mdm2 feedback loop in the presence of DNA damage. The protein p53 activates Mdm2 
and Mdm2 suppresses p53. Morover. the ATM activates the p53 that in turn regulates ATM. 
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Figure 5 - Oscillations of the p53 in the presence of DNA damage 
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showing the oscillations of the p53 in the presence of DNA damage. The bold curve indicates the mean field 
solution and the thin one are results from the Gillespie's simulation. Parameters: k — (0.99, 1, 0.44, 0.69, 0.85, 1.5, 1.5, 0.1, 1.0) h 
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Figure 6 - Power spectrum of the fluctuations of the p53 in the presence of DNA damage 
1 I I I I I 
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frequency f (1/h) 

The points represent the power spectrum obtained after averaging over 200 realizations of the Gillespie's 
algorithm. The continue line is the prediction from the LNA. Note the existence of a characteristic frequency 
of oscillations and a divergence close to zero frequency. Parameters: k = (0.99, 1, 0.44, 0.69, 0.85, 1.5, 1.5, 0.1, 1.0)/i _1 . 
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Figure 7 - p53-Mdm2 model for oncogenic response 




Figure 8 - Oscillations of the p53 after oncogenic response 




30 50 70 90 

time (hours) 

Typical run showing the oscillations of the p53 in the presence of oncogenic signals. The bold curve 
indicates the mean field solution and the thin one are results from the Gillespie's simulation. Parameters: 
k = (1.0, 1.0, 0.44, 0.69, 0.85, 0.5, 0.5, 0.4) 



18 



Figure 9 - Power spectrum of the fluctuations of the p53 after oncogenic response 
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The points represent the power spectrum obtained after averaging over 200 realizations of the Gillespie's 
algorithm. The continue line is the prediction from the LNA. Note the existence of a characteristic frequency 
of oscillations and a divergency close to zero frequency. Parameters: k = (1.0, 1.0, 0.44, 0.69, 0.85, 0.5, 0.5, 0.4) 
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